home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / C / Applications / POV-Ray 3.0.2 / src / SOURCE / HFIELD.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-19  |  44.0 KB  |  2,294 lines  |  [TEXT/CWIE]
  1. /****************************************************************************
  2. *                   hfield.c
  3. *
  4. *    This file implements the height field shape primitive.  The shape is
  5. *    implemented as a collection of triangles which are calculated as
  6. *    needed.  The basic intersection routine first computes the rays
  7. *    intersection with the box marking the limits of the shape, then
  8. *    follows the line from one intersection point to the other, testing the
  9. *    two triangles which form the pixel for an intersection with the ray at
  10. *    each step.
  11. *
  12. *    height field added by Doug Muir with lots of advice and support
  13. *    from David Buck and Drew Wells.
  14. *
  15. *  from Persistence of Vision(tm) Ray Tracer
  16. *  Copyright 1996 Persistence of Vision Team
  17. *---------------------------------------------------------------------------
  18. *  NOTICE: This source code file is provided so that users may experiment
  19. *  with enhancements to POV-Ray and to port the software to platforms other
  20. *  than those supported by the POV-Ray Team.  There are strict rules under
  21. *  which you are permitted to use this file.  The rules are in the file
  22. *  named POVLEGAL.DOC which should be distributed with this file. If
  23. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  24. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  25. *  Forum.  The latest version of POV-Ray may be found there as well.
  26. *
  27. * This program is based on the popular DKB raytracer version 2.12.
  28. * DKBTrace was originally written by David K. Buck.
  29. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  30. *
  31. *****************************************************************************/
  32.  
  33. /****************************************************************************
  34. *
  35. *  Explanation:
  36. *
  37. *    -
  38. *
  39. *  Syntax:
  40. *
  41. *  ---
  42. *
  43. *  Aug 1994 : Merged functions for CSG height fields into functions for
  44. *             non-CSG height fiels. Moved all height field map related
  45. *             data into one data structure. Fixed memory problems. [DB]
  46. *
  47. *  Feb 1995 : Major rewrite of the height field intersection tests. [DB]
  48. *
  49. *****************************************************************************/
  50.  
  51. #include "frame.h"
  52. #include "povray.h"
  53. #include "vector.h"
  54. #include "povproto.h"
  55. #include "hfield.h"
  56. #include "matrices.h"
  57. #include "objects.h"
  58.  
  59.  
  60.  
  61. /*****************************************************************************
  62. * Local preprocessor defines
  63. ******************************************************************************/
  64.  
  65. #define sign(x) (((x) >= 0.0) ? 1.0 : -1.0)
  66.  
  67. #define Get_Height(x, z, HField) ((DBL)(HField)->Data->Map[(z)][(x)])
  68.  
  69. /* Small offest. */
  70.  
  71. #define HFIELD_OFFSET 0.001
  72.  
  73.  
  74.  
  75. /*****************************************************************************
  76. * Static functions
  77. ******************************************************************************/
  78.  
  79. static DBL normalize PARAMS((VECTOR A, VECTOR B));
  80.  
  81. static DBL stretch PARAMS((DBL x));
  82.  
  83. static void smooth_height_field PARAMS((HFIELD *HField, int xsize, int zsize));
  84.  
  85. static int intersect_pixel PARAMS((int x,int z,RAY *Ray,HFIELD *HField,DBL height1,DBL height2));
  86.  
  87. static int add_single_normal PARAMS((HF_VAL **data, int xsize, int zsize,
  88.   int x0, int z0,int x1, int z1,int x2, int z2,VECTOR N));
  89.  
  90. static int  All_HField_Intersections PARAMS((OBJECT *Object,RAY *Ray,ISTACK *Depth_Stack));
  91. static int  Inside_HField PARAMS((VECTOR IPoint,OBJECT *Object));
  92. static void HField_Normal PARAMS((VECTOR Result,OBJECT *Object,INTERSECTION *Inter));
  93. static void Translate_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
  94. static void Rotate_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
  95. static void Scale_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
  96. static void Invert_HField PARAMS((OBJECT *Object));
  97. static void Transform_HField PARAMS((OBJECT *Object,TRANSFORM *Trans));
  98. static void *Copy_HField PARAMS((OBJECT *Object));
  99. static void Destroy_HField PARAMS((OBJECT *Object));
  100.  
  101. static int dda_traversal PARAMS((RAY *Ray, HFIELD *HField, VECTOR Start, HFIELD_BLOCK *Block));
  102. static int block_traversal PARAMS((RAY *Ray, HFIELD *HField, VECTOR Start));
  103. static void build_hfield_blocks PARAMS((HFIELD *HField));
  104.  
  105.  
  106.  
  107. /*****************************************************************************
  108. * Local variables
  109. ******************************************************************************/
  110.  
  111. METHODS HField_Methods =
  112. {
  113.   All_HField_Intersections,
  114.   Inside_HField, HField_Normal,
  115.   Copy_HField, Translate_HField, Rotate_HField,
  116.   Scale_HField, Transform_HField, Invert_HField, Destroy_HField
  117. };
  118.  
  119. static ISTACK *HField_Stack;
  120. static RAY *RRay;
  121. static DBL mindist, maxdist;
  122.  
  123.  
  124. /*****************************************************************************
  125. *
  126. * FUNCTION
  127. *
  128. *   All_HField_Intersections
  129. *
  130. * INPUT
  131. *   
  132. * OUTPUT
  133. *   
  134. * RETURNS
  135. *   
  136. * AUTHOR
  137. *
  138. *   Doug Muir, David Buck, Drew Wells
  139. *   
  140. * DESCRIPTION
  141. *
  142. *   -
  143. *
  144. * CHANGES
  145. *
  146. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  147. *
  148. ******************************************************************************/
  149.  
  150. static int All_HField_Intersections(Object, Ray, Depth_Stack)
  151. OBJECT *Object;
  152. RAY *Ray;
  153. ISTACK *Depth_Stack;
  154. {
  155.   int Side1, Side2;
  156.   VECTOR Start;
  157.   RAY Temp_Ray;
  158.   DBL depth1, depth2;
  159.   HFIELD *HField = (HFIELD *) Object;
  160.  
  161.   Increase_Counter(stats[Ray_HField_Tests]);
  162.  
  163.   MInvTransPoint(Temp_Ray.Initial, Ray->Initial, HField->Trans);
  164.   MInvTransDirection(Temp_Ray.Direction, Ray->Direction, HField->Trans);
  165.  
  166. #ifdef HFIELD_EXTRA_STATS
  167.   Increase_Counter(stats[Ray_HField_Box_Tests]);
  168. #endif
  169.  
  170.   if (!Intersect_Box(&Temp_Ray,HField->bounding_box,&depth1,&depth2,&Side1,&Side2))
  171.   {
  172.     return(FALSE);
  173.   }
  174.  
  175. #ifdef HFIELD_EXTRA_STATS
  176.   Increase_Counter(stats[Ray_HField_Box_Tests_Succeeded]);
  177. #endif
  178.  
  179.   HField->Data->cache_pos = 0;
  180.  
  181.   if (depth1 < EPSILON)
  182.   {
  183.     depth1 = EPSILON;
  184.  
  185.     if (depth1 > depth2)
  186.     {
  187.       return(FALSE);
  188.     }
  189.   }
  190.  
  191.   VEvaluateRay(Start, Temp_Ray.Initial, depth1, Temp_Ray.Direction);
  192.  
  193.   mindist = depth1;
  194.   maxdist = depth2;
  195.  
  196.   HField_Stack = Depth_Stack;
  197.  
  198.   RRay = Ray;
  199.  
  200.   if (block_traversal(&Temp_Ray, HField, Start))
  201.   {
  202.     Increase_Counter(stats[Ray_HField_Tests_Succeeded]);
  203.  
  204.     return(TRUE);
  205.   }
  206.   else
  207.   {
  208.     return(FALSE);
  209.   }
  210. }
  211.  
  212.  
  213.  
  214. /*****************************************************************************
  215. *
  216. * FUNCTION
  217. *
  218. *   Inside_HField
  219. *
  220. * INPUT
  221. *   
  222. * OUTPUT
  223. *   
  224. * RETURNS
  225. *   
  226. * AUTHOR
  227. *
  228. *   Doug Muir, David Buck, Drew Wells
  229. *   
  230. * DESCRIPTION
  231. *
  232. *   -
  233. *
  234. * CHANGES
  235. *
  236. *   -
  237. *
  238. ******************************************************************************/
  239.  
  240. static int Inside_HField (IPoint, Object)
  241. VECTOR IPoint;
  242. OBJECT *Object;
  243. {
  244.   HFIELD *HField = (HFIELD *) Object;
  245.   int px, pz;
  246.   DBL x,z,y1,y2,y3,water, dot1, dot2;
  247.   VECTOR Local_Origin, H_Normal, Test;
  248.  
  249.   MInvTransPoint(Test, IPoint, HField->Trans);
  250.  
  251.   water = HField->bounding_box->bounds[0][Y];
  252.  
  253.   if ((Test[X] < 0.0) || (Test[X] >= HField->bounding_box->bounds[1][X]) ||
  254.       (Test[Z] < 0.0) || (Test[Z] >= HField->bounding_box->bounds[1][Z]))
  255.   {
  256.     return(Test_Flag(HField, INVERTED_FLAG));
  257.   }
  258.  
  259.   if (Test[Y] >= HField->bounding_box->bounds[1][Y])
  260.   {
  261.     return(Test_Flag(HField, INVERTED_FLAG));
  262.   }
  263.  
  264.   if (Test[Y] < water)
  265.   {
  266.     return(!Test_Flag(HField, INVERTED_FLAG));
  267.   }
  268.  
  269.   px = (int)Test[X];
  270.   pz = (int)Test[Z];
  271.  
  272.   x = Test[X] - (DBL)px;
  273.   z = Test[Z] - (DBL)pz;
  274.  
  275.   if ((x+z) < 1.0)
  276.   {
  277.     y1 = max(Get_Height(px,   pz,   HField), water);
  278.     y2 = max(Get_Height(px+1, pz,   HField), water);
  279.     y3 = max(Get_Height(px,   pz+1, HField), water);
  280.  
  281.     Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);
  282.  
  283.     Make_Vector(H_Normal, y1-y2, 1.0, y1-y3);
  284.   }
  285.   else
  286.   {
  287.     px = (int)ceil(Test[X]);
  288.     pz = (int)ceil(Test[Z]);
  289.  
  290.     y1 = max(Get_Height(px,   pz,   HField), water);
  291.     y2 = max(Get_Height(px-1, pz,   HField), water);
  292.     y3 = max(Get_Height(px,   pz-1, HField), water);
  293.  
  294.     Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);
  295.  
  296.     Make_Vector(H_Normal, y2-y1, 1.0, y3-y1);
  297.   }
  298.  
  299.   VDot(dot1, Test, H_Normal);
  300.   VDot(dot2, Local_Origin, H_Normal);
  301.  
  302.   if (dot1 < dot2)
  303.   {
  304.     return(!Test_Flag(HField, INVERTED_FLAG));
  305.   }
  306.  
  307.   return(Test_Flag(HField, INVERTED_FLAG));
  308. }
  309.  
  310.  
  311.  
  312. /*****************************************************************************
  313. *
  314. * FUNCTION
  315. *
  316. *   HField_Normal
  317. *
  318. * INPUT
  319. *   
  320. * OUTPUT
  321. *   
  322. * RETURNS
  323. *   
  324. * AUTHOR
  325. *
  326. *   Doug Muir, David Buck, Drew Wells
  327. *   
  328. * DESCRIPTION
  329. *
  330. *   -
  331. *
  332. * CHANGES
  333. *
  334. *   -
  335. *
  336. ******************************************************************************/
  337.  
  338. static void HField_Normal (Result, Object, Inter)
  339. OBJECT *Object;
  340. VECTOR Result;
  341. INTERSECTION *Inter;
  342. {
  343.   HFIELD *HField = (HFIELD *) Object;
  344.   int px,pz, i;
  345.   DBL x,z,y1,y2,y3,u,v;
  346.   VECTOR Local_Origin;
  347.   VECTOR n[5];
  348.  
  349.   MInvTransPoint(Local_Origin, Inter->IPoint, HField->Trans);
  350.  
  351.   for (i = 0; i < HField->Data->cache_pos; i++)
  352.   {
  353.     if ((Local_Origin[X] == HField->Data->Normal_Cache[i].fx) &&
  354.        (Local_Origin[Z] == HField->Data->Normal_Cache[i].fz))
  355.     {
  356.       Assign_Vector(Result,HField->Data->Normal_Cache[i].normal);
  357.  
  358.       MTransNormal(Result,Result,HField->Trans);
  359.  
  360.       VNormalize(Result,Result);
  361.  
  362.       return;
  363.     }
  364.   }
  365.  
  366.   px = (int)Local_Origin[X];
  367.   pz = (int)Local_Origin[Z];
  368.  
  369.   x = Local_Origin[X] - (DBL)px;
  370.   z = Local_Origin[Z] - (DBL)pz;
  371.  
  372.   if (Test_Flag(HField, SMOOTHED_FLAG))
  373.   {
  374.     n[0][X] = HField->Data->Normals[pz][px][0];
  375.     n[0][Y] = HField->Data->Normals[pz][px][1];
  376.     n[0][Z] = HField->Data->Normals[pz][px][2];
  377.     n[1][X] = HField->Data->Normals[pz][px+1][0];
  378.     n[1][Y] = HField->Data->Normals[pz][px+1][1];
  379.     n[1][Z] = HField->Data->Normals[pz][px+1][2];
  380.     n[2][X] = HField->Data->Normals[pz+1][px][0];
  381.     n[2][Y] = HField->Data->Normals[pz+1][px][1];
  382.     n[2][Z] = HField->Data->Normals[pz+1][px][2];
  383.     n[3][X] = HField->Data->Normals[pz+1][px+1][0];
  384.     n[3][Y] = HField->Data->Normals[pz+1][px+1][1];
  385.     n[3][Z] = HField->Data->Normals[pz+1][px+1][2];
  386.  
  387.     x = stretch(x);
  388.     z = stretch(z);
  389.  
  390.     u = (1.0 - x);
  391.     v = (1.0 - z);
  392.  
  393.     Result[X] = v*(u*n[0][X] + x*n[1][X]) + z*(u*n[2][X] + x*n[3][X]);
  394.     Result[Y] = v*(u*n[0][Y] + x*n[1][Y]) + z*(u*n[2][Y] + x*n[3][Y]);
  395.     Result[Z] = v*(u*n[0][Z] + x*n[1][Z]) + z*(u*n[2][Z] + x*n[3][Z]);
  396.   }
  397.   else
  398.   {
  399.     if ((x+z) <= 1.0)
  400.     {
  401.       /* Lower triangle. */
  402.  
  403.       y1 = Get_Height(px,   pz,   HField);
  404.       y2 = Get_Height(px+1, pz,   HField);
  405.       y3 = Get_Height(px,   pz+1, HField);
  406.  
  407.       Make_Vector(Result, y1-y2, 1.0, y1-y3);
  408.     }
  409.     else
  410.     {
  411.       /* Upper triangle. */
  412.  
  413.       y1 = Get_Height(px+1, pz+1, HField);
  414.       y2 = Get_Height(px,   pz+1, HField);
  415.       y3 = Get_Height(px+1, pz,   HField);
  416.  
  417.       Make_Vector(Result, y2-y1, 1.0, y3-y1);
  418.     }
  419.   }
  420.  
  421.   MTransNormal(Result, Result, HField->Trans);
  422.  
  423.   VNormalize(Result, Result);
  424. }
  425.  
  426.  
  427.  
  428. /*****************************************************************************
  429. *
  430. * FUNCTION
  431. *
  432. *   stretch
  433. *
  434. * INPUT
  435. *   
  436. * OUTPUT
  437. *   
  438. * RETURNS
  439. *   
  440. * AUTHOR
  441. *
  442. *   Doug Muir, David Buck, Drew Wells
  443. *   
  444. * DESCRIPTION
  445. *
  446. *   -
  447. *
  448. * CHANGES
  449. *
  450. *   -
  451. *
  452. ******************************************************************************/
  453.  
  454. static DBL stretch (x)
  455. DBL x;
  456. {
  457.   if (x <= 0.5)
  458.   {
  459.     x = 2.0 * x * x;
  460.   }
  461.   else
  462.   {
  463.     x = 1.0 - (2.0 * (1.0 - x) * (1.0 - x));
  464.   }
  465.  
  466.   return(x);
  467. }
  468.  
  469.  
  470.  
  471. /*****************************************************************************
  472. *
  473. * FUNCTION
  474. *
  475. *   normalize
  476. *
  477. * INPUT
  478. *   
  479. * OUTPUT
  480. *   
  481. * RETURNS
  482. *   
  483. * AUTHOR
  484. *
  485. *   Doug Muir, David Buck, Drew Wells
  486. *   
  487. * DESCRIPTION
  488. *
  489. *   -
  490. *
  491. * CHANGES
  492. *
  493. *   -
  494. *
  495. ******************************************************************************/
  496.  
  497. static DBL normalize(A, B)
  498. VECTOR A, B;
  499. {
  500.   VLength(VTemp, B);
  501.  
  502.   if (fabs(VTemp) > EPSILON)
  503.   {
  504.     VInverseScale(A, B, VTemp);
  505.   }
  506.   else
  507.   {
  508.     Make_Vector(A, 0.0, 1.0, 0.0);
  509.   }
  510.  
  511.   return(VTemp);
  512. }
  513.  
  514.  
  515.  
  516. /*****************************************************************************
  517. *
  518. * FUNCTION
  519. *
  520. *   intersect_pixel
  521. *
  522. * INPUT
  523. *   
  524. * OUTPUT
  525. *   
  526. * RETURNS
  527. *   
  528. * AUTHOR
  529. *
  530. *   Doug Muir, David Buck, Drew Wells
  531. *   
  532. * DESCRIPTION
  533. *
  534. *   -
  535. *
  536. * CHANGES
  537. *
  538. *   -
  539. *
  540. ******************************************************************************/
  541.  
  542. static int intersect_pixel(x, z, Ray, HField, height1, height2)
  543. int x;
  544. int z;
  545. RAY *Ray;
  546. HFIELD *HField;
  547. DBL height1, height2;
  548. {
  549.   int Found;
  550.   DBL dot, depth1, depth2;
  551.   DBL s, t, y1, y2, y3, y4;
  552.   DBL min_y2_y3, max_y2_y3;
  553.   DBL max_height, min_height;
  554.   VECTOR P, N, V1;
  555.  
  556. #ifdef HFIELD_EXTRA_STATS
  557.   Increase_Counter(stats[Ray_HField_Cell_Tests]);
  558. #endif
  559.  
  560.   y1 = Get_Height(x,   z,   HField);
  561.   y2 = Get_Height(x+1, z,   HField);
  562.   y3 = Get_Height(x,   z+1, HField);
  563.   y4 = Get_Height(x+1, z+1, HField);
  564.  
  565.   /* Do we hit this cell at all? */
  566.  
  567.   if (y2 < y3)
  568.   {
  569.     min_y2_y3 = y2;
  570.     max_y2_y3 = y3;
  571.   }
  572.   else
  573.   {
  574.     min_y2_y3 = y3;
  575.     max_y2_y3 = y2;
  576.   }
  577.  
  578.   min_height = min(min(y1, y4), min_y2_y3);
  579.   max_height = max(max(y1, y4), max_y2_y3);
  580.  
  581.   if ((max_height < height1) || (min_height > height2))
  582.   {
  583.     return(FALSE);
  584.   }
  585.  
  586. #ifdef HFIELD_EXTRA_STATS
  587.   Increase_Counter(stats[Ray_HField_Cell_Tests_Succeeded]);
  588. #endif
  589.  
  590.   Found = FALSE;
  591.  
  592.   /* Check if we'll hit first triangle. */
  593.  
  594.   min_height = min(y1, min_y2_y3);
  595.   max_height = max(y1, max_y2_y3);
  596.  
  597.   if ((max_height >= height1) && (min_height <= height2))
  598.   {
  599. #ifdef HFIELD_EXTRA_STATS
  600.     Increase_Counter(stats[Ray_HField_Triangle_Tests]);
  601. #endif
  602.  
  603.     /* Set up triangle. */
  604.  
  605.     Make_Vector(P, (DBL)x, y1, (DBL)z);
  606.  
  607.     /*
  608.      * Calculate the normal vector from:
  609.      *
  610.      * N = V2 x V1, with V1 = <1, y2-y1, 0>, V2 = <0, y3-y1, 1>.
  611.      */
  612.  
  613.     Make_Vector(N, y1-y2, 1.0, y1-y3);
  614.  
  615.     /* Now intersect the triangle. */
  616.  
  617.     VDot(dot, N, Ray->Direction);
  618.  
  619.     if ((dot > EPSILON) || (dot < -EPSILON))
  620.     {
  621.       VSub(V1, P, Ray->Initial);
  622.  
  623.       VDot(depth1, N, V1);
  624.  
  625.       depth1 /= dot;
  626.  
  627.       if ((depth1 >= mindist) && (depth1 <= maxdist))
  628.       {
  629.         s = Ray->Initial[X] + depth1 * Ray->Direction[X] - (DBL)x;
  630.         t = Ray->Initial[Z] + depth1 * Ray->Direction[Z] - (DBL)z;
  631.  
  632.         if ((s >= -0.0001) && (t >= -0.0001) && ((s+t) <= 1.0001))
  633.         {
  634. #ifdef HFIELD_EXTRA_STATS
  635.           Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
  636. #endif
  637.  
  638.           VEvaluateRay(P, RRay->Initial, depth1, RRay->Direction);
  639.  
  640.           if (Point_In_Clip(P, HField->Clip))
  641.           {
  642.             push_entry(depth1, P, (OBJECT *)HField, HField_Stack);
  643.  
  644.             Found = TRUE;
  645.  
  646.             /* Cache normal. */
  647.  
  648.             if (!Test_Flag(HField, SMOOTHED_FLAG))
  649.             {
  650.               if (HField->Data->cache_pos < HF_CACHE_SIZE)
  651.               {
  652.                 Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);
  653.  
  654.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
  655.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;
  656.  
  657.                 HField->Data->cache_pos++;
  658.               }
  659.             }
  660.           }
  661.         }
  662.       }
  663.     }
  664.   }
  665.  
  666.   /* Check if we'll hit second triangle. */
  667.  
  668.   min_height = min(y4, min_y2_y3);
  669.   max_height = max(y4, max_y2_y3);
  670.  
  671.   if ((max_height >= height1) && (min_height <= height2))
  672.   {
  673. #ifdef HFIELD_EXTRA_STATS
  674.     Increase_Counter(stats[Ray_HField_Triangle_Tests]);
  675. #endif
  676.  
  677.     /* Set up triangle. */
  678.  
  679.     Make_Vector(P, (DBL)(x+1), y4, (DBL)(z+1));
  680.  
  681.     /*
  682.      * Calculate the normal vector from:
  683.      *
  684.      * N = V2 x V1, with V1 = <-1, y3-y4, 0>, V2 = <0, y2-y4, -1>.
  685.      */
  686.  
  687.     Make_Vector(N, y3-y4, 1.0, y2-y4);
  688.  
  689.     /* Now intersect the triangle. */
  690.  
  691.     VDot(dot, N, Ray->Direction);
  692.  
  693.     if ((dot > EPSILON) || (dot < -EPSILON))
  694.     {
  695.       VSub(V1, P, Ray->Initial);
  696.  
  697.       VDot(depth2, N, V1);
  698.  
  699.       depth2 /= dot;
  700.  
  701.       if ((depth2 >= mindist) && (depth2 <= maxdist))
  702.       {
  703.         s = Ray->Initial[X] + depth2 * Ray->Direction[X] - (DBL)x;
  704.         t = Ray->Initial[Z] + depth2 * Ray->Direction[Z] - (DBL)z;
  705.  
  706.         if ((s <= 1.0001) && (t <= 1.0001) && ((s+t) >= 0.9999))
  707.         {
  708. #ifdef HFIELD_EXTRA_STATS
  709.           Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
  710. #endif
  711.  
  712.           VEvaluateRay(P, RRay->Initial, depth2, RRay->Direction);
  713.  
  714.           if (Point_In_Clip(P, HField->Clip))
  715.           {
  716.             push_entry(depth2, P, (OBJECT *)HField, HField_Stack);
  717.  
  718.             Found = TRUE;
  719.  
  720.             /* Cache normal. */
  721.  
  722.             if (!Test_Flag(HField, SMOOTHED_FLAG))
  723.             {
  724.               if (HField->Data->cache_pos < HF_CACHE_SIZE)
  725.               {
  726.                 Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);
  727.  
  728.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
  729.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;
  730.  
  731.                 HField->Data->cache_pos++;
  732.               }
  733.             }
  734.           }
  735.         }
  736.       }
  737.     }
  738.   }
  739.  
  740.   return(Found);
  741. }
  742.  
  743.  
  744.  
  745. /*****************************************************************************
  746. *
  747. * FUNCTION
  748. *
  749. *   add_single_normal
  750. *
  751. * INPUT
  752. *   
  753. * OUTPUT
  754. *   
  755. * RETURNS
  756. *   
  757. * AUTHOR
  758. *
  759. *   Doug Muir, David Buck, Drew Wells
  760. *   
  761. * DESCRIPTION
  762. *
  763. *   -
  764. *
  765. * CHANGES
  766. *
  767. *   -
  768. *
  769. ******************************************************************************/
  770.  
  771. static int add_single_normal(data, xsize, zsize, x0, z0, x1, z1, x2, z2, N)
  772. HF_VAL **data;
  773. int xsize,zsize,x0,z0,x1,z1,x2,z2;
  774. VECTOR N;
  775. {
  776.   VECTOR v0, v1, v2, t0, t1, Nt;
  777.  
  778.   if ((x0 < 0) || (z0 < 0) ||
  779.       (x1 < 0) || (z1 < 0) ||
  780.       (x2 < 0) || (z2 < 0) ||
  781.       (x0 > xsize) || (z0 > zsize) ||
  782.       (x1 > xsize) || (z1 > zsize) ||
  783.       (x2 > xsize) || (z2 > zsize))
  784.   {
  785.     return(0);
  786.   }
  787.   else
  788.   {
  789.     Make_Vector(v0, x0, (DBL)data[z0][x0], z0);
  790.     Make_Vector(v1, x1, (DBL)data[z1][x1], z1);
  791.     Make_Vector(v2, x2, (DBL)data[z2][x2], z2);
  792.  
  793.     VSub(t0, v2, v0);
  794.     VSub(t1, v1, v0);
  795.  
  796.     VCross(Nt, t0, t1);
  797.  
  798.     normalize(Nt, Nt);
  799.  
  800.     if (Nt[Y] < 0.0)
  801.     {
  802.       VScaleEq(Nt, -1.0);
  803.     }
  804.  
  805.     VAddEq(N, Nt);
  806.  
  807.     return(1);
  808.   }
  809. }
  810.  
  811.  
  812.  
  813. /*****************************************************************************
  814. *
  815. * FUNCTION
  816. *
  817. *   smooth_height_field
  818. *
  819. * INPUT
  820. *   
  821. * OUTPUT
  822. *   
  823. * RETURNS
  824. *   
  825. * AUTHOR
  826. *
  827. *   Doug Muir, David Buck, Drew Wells
  828. *   
  829. * DESCRIPTION
  830. *
  831. *   Given a height field that only contains an elevation grid, this
  832. *   routine will walk through the data and produce averaged normals
  833. *   for all points on the grid.
  834. *
  835. * CHANGES
  836. *
  837. *   -
  838. *
  839. ******************************************************************************/
  840.  
  841. static void smooth_height_field(HField, xsize, zsize)
  842. HFIELD *HField;
  843. int xsize;
  844. int zsize;
  845. {
  846.   int i, j, k;
  847.   VECTOR N;
  848.   HF_VAL **map = HField->Data->Map;
  849.  
  850.   /* First off, allocate all the memory needed to store the normal information */
  851.  
  852.   HField->Data->Normals_Height = zsize+1;
  853.  
  854.   HField->Data->Normals = (HF_Normals **)POV_MALLOC((zsize+1)*sizeof(HF_Normals *), "height field normals");
  855.  
  856.   for (i = 0; i <= zsize; i++)
  857.   {
  858.     HField->Data->Normals[i] = (HF_Normals *)POV_MALLOC((xsize+1)*sizeof(HF_Normals), "height field normals");
  859.   }
  860.  
  861.   /*
  862.    * For now we will do it the hard way - by generating the normals
  863.    * individually for each elevation point.
  864.    */
  865.  
  866.   for (i = 0; i < zsize; i++)
  867.   {
  868.     COOPERATE_1
  869.  
  870.     for (j = 0; j < xsize; j++)
  871.     {
  872.       Make_Vector(N, 0.0, 0.0, 0.0);
  873.  
  874.       k = 0;
  875.  
  876.       k += add_single_normal(map, xsize, zsize, j, i, j+1, i, j, i+1, N);
  877.       k += add_single_normal(map, xsize, zsize, j, i, j, i+1, j-1, i, N);
  878.       k += add_single_normal(map, xsize, zsize, j, i, j-1, i, j, i-1, N);
  879.       k += add_single_normal(map, xsize, zsize, j, i, j, i-1, j+1, i, N);
  880.  
  881.       if (k == 0)
  882.       {
  883.         Error ("Failed to find any normals at: (%d, %d).\n", i, j);
  884.       }
  885.  
  886.       normalize(N, N);
  887.  
  888.       HField->Data->Normals[i][j][0] = (short)(32767 * N[X]);
  889.       HField->Data->Normals[i][j][1] = (short)(32767 * N[Y]);
  890.       HField->Data->Normals[i][j][2] = (short)(32767 * N[Z]);
  891.     }
  892.   }
  893. }
  894.  
  895.  
  896.  
  897. /*****************************************************************************
  898. *
  899. * FUNCTION
  900. *
  901. *   Compute_HField
  902. *
  903. * INPUT
  904. *
  905. * OUTPUT
  906. *
  907. * RETURNS
  908. *
  909. * AUTHOR
  910. *
  911. *   Doug Muir, David Buck, Drew Wells
  912. *
  913. * DESCRIPTION
  914. *
  915. *   Copy image data into height field map. Create bounding blocks
  916. *   for the block traversal. Calculate normals for smoothed height fields.
  917. *
  918. * CHANGES
  919. *
  920. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  921. *
  922. ******************************************************************************/
  923.  
  924. void Compute_HField(HField, Image)
  925. HFIELD *HField;
  926. IMAGE *Image;
  927. {
  928.   int x, z, max_x, max_z;
  929.   int temp1 = 0, temp2 = 0;
  930.   HF_VAL min_y, max_y, temp_y;
  931.  
  932.   /* Get height field map size. */
  933.  
  934.   max_x = Image->iwidth;
  935.  
  936.   if (Image->File_Type == POT_FILE)
  937.   {
  938.     max_x = max_x / 2;
  939.   }
  940.  
  941.   max_z = Image->iheight;
  942.  
  943.   /* Allocate memory for map. */
  944.  
  945.   HField->Data->Map = (HF_VAL **)POV_MALLOC(max_z*sizeof(HF_VAL *), "height field");
  946.  
  947.   for (z = 0; z < max_z; z++)
  948.   {
  949.     HField->Data->Map[z] = (HF_VAL *)POV_MALLOC(max_x*sizeof(HF_VAL), "height field");
  950.   }
  951.  
  952.   /* Copy map. */
  953.  
  954.   min_y = 65535L;
  955.   max_y = 0;
  956.  
  957.   for (z = 0; z < max_z; z++)
  958.   {
  959.     COOPERATE_1
  960.     for (x = 0; x < max_x; x++)
  961.     {
  962.  
  963.       switch (Image->File_Type)
  964.       {
  965.         case GIF_FILE:
  966.  
  967.           temp1 = Image->data.map_lines[max_z - z - 1][x];
  968.           temp2 = 0;
  969.  
  970.           break;
  971.  
  972.         case POT_FILE:
  973.  
  974.           temp1 = Image->data.map_lines[max_z - z - 1][x];
  975.           temp2 = Image->data.map_lines[max_z - z - 1][x + max_x];
  976.  
  977.           break;
  978.  
  979.  
  980.         case PPM_FILE:
  981.  
  982.           temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
  983.           temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];
  984.  
  985.           break;
  986.  
  987.         case PGM_FILE:
  988.         case TGA_FILE:
  989.         case PNG_FILE:
  990.  
  991.           if (Image->Colour_Map == NULL)
  992.           {
  993.             temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
  994.             temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];
  995.           }
  996.           else
  997.           {
  998.             temp1 = Image->data.map_lines[max_z - z - 1][x];
  999.             temp2 = 0;
  1000.           }
  1001.  
  1002.           break;
  1003.  
  1004.         default:
  1005.  
  1006.           Error("Unknown image type in Compute_HField().\n");
  1007.       }
  1008.  
  1009.       temp_y = (HF_VAL)(256*temp1 + temp2);
  1010.  
  1011.       HField->Data->Map[z][x] = temp_y;
  1012.  
  1013.       min_y = min(min_y, temp_y);
  1014.       max_y = max(max_y, temp_y);
  1015.     }
  1016.   }
  1017.  
  1018.   /* Resize bounding box. */
  1019.  
  1020.   HField->Data->min_y = min_y;
  1021.   HField->Data->max_y = max_y;
  1022.  
  1023.   HField->bounding_box->bounds[0][Y] = max((DBL)min_y, HField->bounding_box->bounds[0][Y]) - HFIELD_OFFSET;
  1024.   HField->bounding_box->bounds[1][Y] = (DBL)max_y + HFIELD_OFFSET;
  1025.  
  1026.   /* Compute smoothed height field. */
  1027.  
  1028.   if (Test_Flag(HField, SMOOTHED_FLAG))
  1029.   {
  1030.     smooth_height_field(HField, max_x-1, max_z-1);
  1031.   }
  1032.  
  1033.   HField->Data->max_x = max_x-2;
  1034.   HField->Data->max_z = max_z-2;
  1035.  
  1036.   build_hfield_blocks(HField);
  1037. }
  1038.  
  1039.  
  1040.  
  1041. /*****************************************************************************
  1042. *
  1043. * FUNCTION
  1044. *
  1045. *   build_hfield_blocks
  1046. *
  1047. * INPUT
  1048. *   
  1049. * OUTPUT
  1050. *   
  1051. * RETURNS
  1052. *   
  1053. * AUTHOR
  1054. *
  1055. *   Dieter Bayer
  1056. *   
  1057. * DESCRIPTION
  1058. *
  1059. *   Create the bounding block hierarchy used by the block traversal.
  1060. *
  1061. * CHANGES
  1062. *
  1063. *   Feb 1995 : Creation.
  1064. *
  1065. ******************************************************************************/
  1066.  
  1067. static void build_hfield_blocks(HField)
  1068. HFIELD *HField;
  1069. {
  1070.   int x, z, nx, nz, wx, wz;
  1071.   int i, j;
  1072.   int xmin, xmax, zmin, zmax;
  1073.   DBL y, ymin, ymax, water;
  1074.  
  1075.   /* Get block size. */
  1076.  
  1077.   nx = max(1, (int)(sqrt((DBL)HField->Data->max_x)));
  1078.   nz = max(1, (int)(sqrt((DBL)HField->Data->max_z)));
  1079.  
  1080.   /* Get dimensions of sub-block. */
  1081.  
  1082.   wx = (int)ceil((DBL)(HField->Data->max_x + 2) / (DBL)nx);
  1083.   wz = (int)ceil((DBL)(HField->Data->max_z + 2) / (DBL)nz);
  1084.  
  1085.   /* Increase number of sub-blocks if necessary. */
  1086.  
  1087.   if (nx * wx < HField->Data->max_x + 2)
  1088.   {
  1089.     nx++;
  1090.   }
  1091.  
  1092.   if (nz * wz < HField->Data->max_z + 2)
  1093.   {
  1094.     nz++;
  1095.   }
  1096.  
  1097.   if (!Test_Flag(HField, HIERARCHY_FLAG) || ((nx == 1) && (nz == 1)))
  1098.   {
  1099.     /* We don't want a bounding hierarchy. Just use one block. */
  1100.  
  1101.     HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(sizeof(HFIELD_BLOCK *), "height field blocks");
  1102.  
  1103.     HField->Data->Block[0] = (HFIELD_BLOCK *)POV_MALLOC(sizeof(HFIELD_BLOCK), "height field blocks");
  1104.  
  1105.     HField->Data->Block[0][0].xmin = 0;
  1106.     HField->Data->Block[0][0].xmax = HField->Data->max_x;
  1107.     HField->Data->Block[0][0].zmin = 0;
  1108.     HField->Data->Block[0][0].zmax = HField->Data->max_z;
  1109.  
  1110.     HField->Data->Block[0][0].ymin = HField->bounding_box->bounds[0][Y];
  1111.     HField->Data->Block[0][0].ymax = HField->bounding_box->bounds[1][Y];
  1112.  
  1113.     HField->Data->block_max_x = 1;
  1114.     HField->Data->block_max_z = 1;
  1115.  
  1116.     HField->Data->block_width_x = HField->Data->max_x + 2;
  1117.     HField->Data->block_width_z = HField->Data->max_y + 2;
  1118.  
  1119. /*
  1120.     Debug_Info("\nHeight field: %d x %d (1 x 1 blocks)", HField->Data->max_x+2, HField->Data->max_z+2);
  1121. */
  1122.  
  1123.     return;
  1124.   }
  1125.  
  1126. /*
  1127.   Debug_Info("\nHeight field: %d x %d (%d x %d blocks)", HField->Data->max_x+2, HField->Data->max_z+2, nx, nz);
  1128. */
  1129.  
  1130.   /* Allocate memory for blocks. */
  1131.  
  1132.   HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(nz*sizeof(HFIELD_BLOCK *), "height field blocks");
  1133.  
  1134.   /* Store block information. */
  1135.  
  1136.   HField->Data->block_max_x = nx;
  1137.   HField->Data->block_max_z = nz;
  1138.  
  1139.   HField->Data->block_width_x = wx;
  1140.   HField->Data->block_width_z = wz;
  1141.  
  1142.   water = HField->bounding_box->bounds[0][Y];
  1143.  
  1144.   for (z = 0; z < nz; z++)
  1145.   {
  1146.     COOPERATE_1
  1147.  
  1148.     HField->Data->Block[z] = (HFIELD_BLOCK *)POV_MALLOC(nx*sizeof(HFIELD_BLOCK), "height field blocks");
  1149.  
  1150.     for (x = 0; x < nx; x++)
  1151.     {
  1152.       /* Get block's borders. */
  1153.  
  1154.       xmin = x * wx;
  1155.       zmin = z * wz;
  1156.  
  1157.       xmax = min((x + 1) * wx - 1, HField->Data->max_x);
  1158.       zmax = min((z + 1) * wz - 1, HField->Data->max_z);
  1159.  
  1160.       /* Find min. and max. height in current block. */
  1161.  
  1162.       ymin = BOUND_HUGE;
  1163.       ymax = -BOUND_HUGE;
  1164.  
  1165.       for (i = xmin; i <= xmax+1; i++)
  1166.       {
  1167.         for (j = zmin; j <= zmax+1; j++)
  1168.         {
  1169.           y = Get_Height(i, j, HField);
  1170.  
  1171.           ymin = min(ymin, y);
  1172.           ymax = max(ymax, y);
  1173.         }
  1174.       }
  1175.  
  1176.       /* Store block's borders. */
  1177.  
  1178.       HField->Data->Block[z][x].xmin = xmin;
  1179.       HField->Data->Block[z][x].xmax = xmax;
  1180.       HField->Data->Block[z][x].zmin = zmin;
  1181.       HField->Data->Block[z][x].zmax = zmax;
  1182.  
  1183.       HField->Data->Block[z][x].ymin = max(ymin, water) - HFIELD_OFFSET;
  1184.       HField->Data->Block[z][x].ymax = ymax + HFIELD_OFFSET;
  1185.     }
  1186.   }
  1187. }
  1188.  
  1189.  
  1190.  
  1191. /*****************************************************************************
  1192. *
  1193. * FUNCTION
  1194. *
  1195. *   Translate_HField
  1196. *
  1197. * INPUT
  1198. *   
  1199. * OUTPUT
  1200. *   
  1201. * RETURNS
  1202. *   
  1203. * AUTHOR
  1204. *
  1205. *   Doug Muir, David Buck, Drew Wells
  1206. *
  1207. * DESCRIPTION
  1208. *
  1209. *   -
  1210. *
  1211. * CHANGES
  1212. *
  1213. *   -
  1214. *
  1215. ******************************************************************************/
  1216.  
  1217. static void Translate_HField (Object, Vector, Trans)
  1218. OBJECT *Object;
  1219. VECTOR Vector;
  1220. TRANSFORM *Trans;
  1221. {
  1222.   Transform_HField(Object, Trans);
  1223. }
  1224.  
  1225.  
  1226.  
  1227. /*****************************************************************************
  1228. *
  1229. * FUNCTION
  1230. *
  1231. *   Rotate_HField
  1232. *
  1233. * INPUT
  1234. *
  1235. * OUTPUT
  1236. *
  1237. * RETURNS
  1238. *
  1239. * AUTHOR
  1240. *
  1241. *   Doug Muir, David Buck, Drew Wells
  1242. *
  1243. * DESCRIPTION
  1244. *
  1245. *   -
  1246. *
  1247. * CHANGES
  1248. *
  1249. *   -
  1250. *
  1251. ******************************************************************************/
  1252.  
  1253. static void Rotate_HField (Object, Vector, Trans)
  1254. OBJECT *Object;
  1255. VECTOR Vector;
  1256. TRANSFORM *Trans;
  1257. {
  1258.   Transform_HField(Object, Trans);
  1259. }
  1260.  
  1261.  
  1262.  
  1263. /*****************************************************************************
  1264. *
  1265. * FUNCTION
  1266. *
  1267. *   Scale_HField
  1268. *
  1269. * INPUT
  1270. *   
  1271. * OUTPUT
  1272. *   
  1273. * RETURNS
  1274. *   
  1275. * AUTHOR
  1276. *
  1277. *   Doug Muir, David Buck, Drew Wells
  1278. *   
  1279. * DESCRIPTION
  1280. *
  1281. *   -
  1282. *
  1283. * CHANGES
  1284. *
  1285. *   -
  1286. *
  1287. ******************************************************************************/
  1288.  
  1289. static void Scale_HField (Object, Vector, Trans)
  1290. OBJECT *Object;
  1291. VECTOR Vector;
  1292. TRANSFORM *Trans;
  1293. {
  1294.   Transform_HField(Object, Trans);
  1295. }
  1296.  
  1297.  
  1298.  
  1299. /*****************************************************************************
  1300. *
  1301. * FUNCTION
  1302. *
  1303. *   Invert_HField
  1304. *
  1305. * INPUT
  1306. *   
  1307. * OUTPUT
  1308. *   
  1309. * RETURNS
  1310. *   
  1311. * AUTHOR
  1312. *
  1313. *   Doug Muir, David Buck, Drew Wells
  1314. *   
  1315. * DESCRIPTION
  1316. *
  1317. *   -
  1318. *
  1319. * CHANGES
  1320. *
  1321. *   -
  1322. *
  1323. ******************************************************************************/
  1324.  
  1325. static void Invert_HField (Object)
  1326. OBJECT *Object;
  1327. {
  1328.   Invert_Flag(Object, INVERTED_FLAG);
  1329. }
  1330.  
  1331.  
  1332.  
  1333. /*****************************************************************************
  1334. *
  1335. * FUNCTION
  1336. *
  1337. *   Transform_HField
  1338. *
  1339. * INPUT
  1340. *   
  1341. * OUTPUT
  1342. *   
  1343. * RETURNS
  1344. *   
  1345. * AUTHOR
  1346. *
  1347. *   Doug Muir, David Buck, Drew Wells
  1348. *   
  1349. * DESCRIPTION
  1350. *
  1351. *   -
  1352. *
  1353. * CHANGES
  1354. *
  1355. *   -
  1356. *
  1357. ******************************************************************************/
  1358.  
  1359. static void Transform_HField (Object, Trans)
  1360. OBJECT *Object;
  1361. TRANSFORM *Trans;
  1362. {
  1363.   Compose_Transforms(((HFIELD *)Object)->Trans, Trans);
  1364.  
  1365.   Compute_HField_BBox((HFIELD *)Object);
  1366. }
  1367.  
  1368.  
  1369.  
  1370. /*****************************************************************************
  1371. *
  1372. * FUNCTION
  1373. *
  1374. *   Create_HField
  1375. *
  1376. * INPUT
  1377. *
  1378. * OUTPUT
  1379. *   
  1380. * RETURNS
  1381. *   
  1382. * AUTHOR
  1383. *
  1384. *   Doug Muir, David Buck, Drew Wells
  1385. *   
  1386. * DESCRIPTION
  1387. *
  1388. *   Allocate and intialize a height field.
  1389. *
  1390. * CHANGES
  1391. *
  1392. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  1393. *
  1394. ******************************************************************************/
  1395.  
  1396. HFIELD *Create_HField()
  1397. {
  1398.   HFIELD *New;
  1399.  
  1400.   /* Allocate height field. */
  1401.  
  1402.   New = (HFIELD *)POV_MALLOC(sizeof(HFIELD), "height field");
  1403.  
  1404.   INIT_OBJECT_FIELDS(New, HFIELD_OBJECT, &HField_Methods)
  1405.  
  1406.   /* Always uses Trans so always create one. */
  1407.  
  1408.   New->Trans = Create_Transform();
  1409.  
  1410.   New->bounding_box = Create_Box();
  1411.  
  1412.   /* Allocate height field data. */
  1413.  
  1414.   New->Data = (HFIELD_DATA *)POV_MALLOC(sizeof(HFIELD_DATA), "height field");
  1415.  
  1416.   New->Data->References = 1;
  1417.  
  1418.   New->Data->cache_pos = 0;
  1419.  
  1420.   New->Data->Normals_Height = 0;
  1421.  
  1422.   New->Data->Map     = NULL;
  1423.   New->Data->Normals = NULL;
  1424.  
  1425.   New->Data->max_x = 0;
  1426.   New->Data->max_z = 0;
  1427.  
  1428.   New->Data->block_max_x = 0;
  1429.   New->Data->block_max_z = 0;
  1430.  
  1431.   New->Data->block_width_x = 0;
  1432.   New->Data->block_width_z = 0;
  1433.  
  1434.   Set_Flag(New, HIERARCHY_FLAG);
  1435.  
  1436.   return(New);
  1437. }
  1438.  
  1439.  
  1440.  
  1441. /*****************************************************************************
  1442. *
  1443. * FUNCTION
  1444. *
  1445. *   Copy_HField
  1446. *
  1447. * INPUT
  1448. *   
  1449. * OUTPUT
  1450. *   
  1451. * RETURNS
  1452. *   
  1453. * AUTHOR
  1454. *
  1455. *   Doug Muir, David Buck, Drew Wells
  1456. *   
  1457. * DESCRIPTION
  1458. *
  1459. *   NOTE: The height field data is not copied, only the number of references
  1460. *         is counted, so that Destray_HField() knows if it can be destroyed.
  1461. *
  1462. * CHANGES
  1463. *
  1464. *   -
  1465. *
  1466. ******************************************************************************/
  1467.  
  1468. static void *Copy_HField(Object)
  1469. OBJECT *Object;
  1470. {
  1471.   HFIELD *New;
  1472.  
  1473.   New = Create_HField();
  1474.  
  1475.   /* Destroy obsolete things created in Create_HField(). */
  1476.  
  1477.   Destroy_Transform(New->Trans);
  1478.  
  1479.   Destroy_Box((OBJECT *)(New->bounding_box));
  1480.  
  1481.   POV_FREE (New->Data);
  1482.  
  1483.   /* Copy height field. */
  1484.  
  1485.   *New = *((HFIELD *)Object);
  1486.  
  1487.   New->Trans = Copy_Transform(((HFIELD *)Object)->Trans);
  1488.  
  1489.   New->bounding_box = Copy_Box((OBJECT *)(((HFIELD *)Object)->bounding_box));
  1490.  
  1491.   New->Data->References++;
  1492.  
  1493.   return(New);
  1494. }
  1495.  
  1496.  
  1497.  
  1498. /*****************************************************************************
  1499. *
  1500. * FUNCTION
  1501. *
  1502. *   Destroy_HField
  1503. *
  1504. * INPUT
  1505. *   
  1506. * OUTPUT
  1507. *   
  1508. * RETURNS
  1509. *   
  1510. * AUTHOR
  1511. *
  1512. *   Doug Muir, David Buck, Drew Wells
  1513. *   
  1514. * DESCRIPTION
  1515. *
  1516. *   NOTE: The height field data is destroyed if it's no longer
  1517. *         used by any copy.
  1518. *
  1519. * CHANGES
  1520. *
  1521. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  1522. *
  1523. ******************************************************************************/
  1524.  
  1525. static void Destroy_HField (Object)
  1526. OBJECT *Object;
  1527. {
  1528.   int i;
  1529.   HFIELD *HField = (HFIELD *)Object;
  1530.  
  1531.   Destroy_Transform(HField->Trans);
  1532.  
  1533.   Destroy_Box((OBJECT *)(HField->bounding_box));
  1534.  
  1535.   if (--(HField->Data->References) == 0)
  1536.   {
  1537.     if (HField->Data->Map != NULL)
  1538.     {
  1539.       for (i = 0; i < HField->Data->max_z+2; i++)
  1540.       {
  1541.         if (HField->Data->Map[i] != NULL)
  1542.         {
  1543.           POV_FREE (HField->Data->Map[i]);
  1544.         }
  1545.       }
  1546.  
  1547.       POV_FREE (HField->Data->Map);
  1548.     }
  1549.  
  1550.     if (HField->Data->Normals != NULL)
  1551.     {
  1552.       for (i = 0; i < HField->Data->Normals_Height; i++)
  1553.       {
  1554.         POV_FREE (HField->Data->Normals[i]);
  1555.       }
  1556.  
  1557.       POV_FREE (HField->Data->Normals);
  1558.     }
  1559.  
  1560.     if (HField->Data->Block != NULL)
  1561.     {
  1562.       for (i = 0; i < HField->Data->block_max_z; i++)
  1563.       {
  1564.         POV_FREE(HField->Data->Block[i]);
  1565.       }
  1566.  
  1567.       POV_FREE(HField->Data->Block);
  1568.     }
  1569.  
  1570.     POV_FREE (HField->Data);
  1571.   }
  1572.  
  1573.   POV_FREE (Object);
  1574. }
  1575.  
  1576.  
  1577.  
  1578. /*****************************************************************************
  1579. *
  1580. * FUNCTION
  1581. *
  1582. *   Compute_HField_BBox
  1583. *
  1584. * INPUT
  1585. *
  1586. *   HField - Height field
  1587. *   
  1588. * OUTPUT
  1589. *
  1590. *   HField
  1591. *   
  1592. * RETURNS
  1593. *   
  1594. * AUTHOR
  1595. *
  1596. *   Dieter Bayer
  1597. *   
  1598. * DESCRIPTION
  1599. *
  1600. *   Calculate the bounding box of a height field.
  1601. *
  1602. * CHANGES
  1603. *
  1604. *   Aug 1994 : Creation.
  1605. *
  1606. ******************************************************************************/
  1607.  
  1608. void Compute_HField_BBox(HField)
  1609. HFIELD *HField;
  1610. {
  1611.   Assign_BBox_Vect(HField->BBox.Lower_Left, HField->bounding_box->bounds[0]);
  1612.  
  1613.   VSub (HField->BBox.Lengths, HField->bounding_box->bounds[1], HField->bounding_box->bounds[0]);
  1614.  
  1615.   if (HField->Trans != NULL)
  1616.   {
  1617.     Recompute_BBox(&HField->BBox, HField->Trans);
  1618.   }
  1619. }
  1620.  
  1621.  
  1622.  
  1623. /*****************************************************************************
  1624. *
  1625. * FUNCTION
  1626. *
  1627. *   dda_traversal
  1628. *
  1629. * INPUT
  1630. *
  1631. *   Ray    - Current ray
  1632. *   HField - Height field
  1633. *   Start  - Start point for the walk
  1634. *   Block  - Sub-block of the height field to traverse
  1635. *   
  1636. * OUTPUT
  1637. *   
  1638. * RETURNS
  1639. *
  1640. *   int - TRUE if intersection was found
  1641. *   
  1642. * AUTHOR
  1643. *
  1644. *   Dieter Bayer
  1645. *   
  1646. * DESCRIPTION
  1647. *
  1648. *   Traverse the grid cell of the height field using a modified DDA.
  1649. *
  1650. *   Based on the following article:
  1651. *
  1652. *     Musgrave, F. Kenton, "Grid Tracing: Fast Ray Tracing for Height
  1653. *     Fields", Research Report YALEU-DCS-RR-39, Yale University, July 1988
  1654. *
  1655. *   You should note that there are (n-1) x (m-1) grid cells in a height
  1656. *   field of (image) size n x m. A grid cell (i,j), 0 <= i <= n-1,
  1657. *   0 <= j <= m-1, extends from x = i to x = i + 1 - epsilon and
  1658. *   y = j to y = j + 1 -epsilon.
  1659. *
  1660. * CHANGES
  1661. *
  1662. *   Feb 1995 : Creation.
  1663. *
  1664. ******************************************************************************/
  1665.  
  1666. static int dda_traversal(Ray, HField, Start, Block)
  1667. RAY *Ray;
  1668. HFIELD *HField;
  1669. VECTOR Start;
  1670. HFIELD_BLOCK *Block;
  1671. {
  1672.   int found;
  1673.   int xmin, xmax, zmin, zmax;
  1674.   int x, z, signx, signz;
  1675.   DBL ymin, ymax, y1, y2;
  1676.   DBL px, pz, dx, dy, dz;
  1677.   DBL delta, error, x0, z0;
  1678.   DBL neary, fary, deltay;
  1679.  
  1680.   /* Setup DDA. */
  1681.  
  1682.   found = FALSE;
  1683.  
  1684.   px = Start[X];
  1685.   pz = Start[Z];
  1686.  
  1687.   /* Get dimensions of current block. */
  1688.  
  1689.   xmin = Block->xmin;
  1690.   xmax = min(Block->xmax + 1, HField->Data->max_x);
  1691.   zmin = Block->zmin;
  1692.   zmax = min(Block->zmax + 1, HField->Data->max_z);
  1693.  
  1694.   ymin = min(Start[Y], Block->ymin) - EPSILON;
  1695.   ymax = max(Start[Y], Block->ymax) + EPSILON;
  1696.  
  1697.   /* Check for illegal grid values (caused by numerical inaccuracies). */
  1698.  
  1699.   if (px < (DBL)xmin)
  1700.   {
  1701.     if (px < (DBL)xmin - HFIELD_OFFSET)
  1702.     {
  1703.       Debug_Info("Illegal grid value in dda_traversal().\n");
  1704.  
  1705.       return(FALSE);
  1706.     }
  1707.     else
  1708.     {
  1709.       px = (DBL)xmin;
  1710.     }
  1711.   }
  1712.   else
  1713.   {
  1714.     if (px > (DBL)xmax + 1.0 - EPSILON)
  1715.     {
  1716.       if (px > (DBL)xmax + 1.0 + EPSILON)
  1717.       {
  1718.         Debug_Info("Illegal grid value in dda_traversal().\n");
  1719.  
  1720.         return(FALSE);
  1721.       }
  1722.       else
  1723.       {
  1724.         px = (DBL)xmax + 1.0 - EPSILON;
  1725.       }
  1726.     }
  1727.   }
  1728.  
  1729.   if (pz < (DBL)zmin)
  1730.   {
  1731.     if (pz < (DBL)zmin - HFIELD_OFFSET)
  1732.     {
  1733.       Debug_Info("Illegal grid value in dda_traversal().\n");
  1734.  
  1735.       return(FALSE);
  1736.     }
  1737.     else
  1738.     {
  1739.       pz = (DBL)zmin;
  1740.     }
  1741.   }
  1742.   else
  1743.   {
  1744.     if (pz > (DBL)zmax + 1.0 - EPSILON)
  1745.     {
  1746.       if (pz > (DBL)zmax + 1.0 + EPSILON)
  1747.       {
  1748.         Debug_Info("Illegal grid value in dda_traversal().\n");
  1749.  
  1750.         return(FALSE);
  1751.       }
  1752.       else
  1753.       {
  1754.         pz = (DBL)zmax + 1.0 - EPSILON;
  1755.       }
  1756.     }
  1757.   }
  1758.  
  1759.   dx = Ray->Direction[X];
  1760.   dy = Ray->Direction[Y];
  1761.   dz = Ray->Direction[Z];
  1762.  
  1763.   /*
  1764.    * Here comes the DDA algorithm.
  1765.    */
  1766.  
  1767.   /* Choose algorithm depending on the driving axis. */
  1768.  
  1769.   if (fabs(dx) >= fabs(dz))
  1770.   {
  1771.     /*
  1772.      * X-axis is driving axis.
  1773.      */
  1774.  
  1775.     delta = fabs(dz / dx);
  1776.  
  1777.     x = (int)px;
  1778.     z = (int)pz;
  1779.  
  1780.     x0 = px - floor(px);
  1781.     z0 = pz - floor(pz);
  1782.  
  1783.     signx = sign(dx);
  1784.     signz = sign(dz);
  1785.  
  1786.     /* Get initial error. */
  1787.  
  1788.     if (dx >= 0.0)
  1789.     {
  1790.       if (dz >= 0.0)
  1791.       {
  1792.         error = z0 + delta * (1.0 - x0) - 1.0;
  1793.       }
  1794.       else
  1795.       {
  1796.         error = -(z0 - delta * (1.0 - x0));
  1797.       }
  1798.     }
  1799.     else
  1800.     {
  1801.       if (dz >= 0.0)
  1802.       {
  1803.         error = z0 + delta * x0 - 1.0;
  1804.       }
  1805.       else
  1806.       {
  1807.         error = -(z0 - delta * x0);
  1808.       }
  1809.     }
  1810.  
  1811.     /* Get y differential. */
  1812.  
  1813.     deltay = dy / fabs(dx);
  1814.  
  1815.     if (dx >= 0.0)
  1816.     {
  1817.       neary = Start[Y] - x0 * deltay;
  1818.  
  1819.       fary = neary + deltay;
  1820.     }
  1821.     else
  1822.     {
  1823.       neary = Start[Y] - (1.0 - x0) * deltay;
  1824.  
  1825.       fary = neary + deltay;
  1826.     }
  1827.  
  1828.     /* Step through the cells. */
  1829.  
  1830.     do
  1831.     {
  1832.       if (neary < fary)
  1833.       {
  1834.         y1 = neary;
  1835.         y2 = fary;
  1836.       }
  1837.       else
  1838.       {
  1839.         y1 = fary;
  1840.         y2 = neary;
  1841.       }
  1842.  
  1843.       if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1844.       {
  1845.         if (HField->Type & IS_CHILD_OBJECT)
  1846.         {
  1847.           found = TRUE;
  1848.         }
  1849.         else
  1850.         {
  1851.           return(TRUE);
  1852.         }
  1853.       }
  1854.  
  1855.       if (error > EPSILON)
  1856.       {
  1857.         z += signz;
  1858.  
  1859.         if ((z < zmin) || (z > zmax))
  1860.         {
  1861.           break;
  1862.         }
  1863.         else
  1864.         {
  1865.           if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1866.           {
  1867.             if (HField->Type & IS_CHILD_OBJECT)
  1868.             {
  1869.               found = TRUE;
  1870.             }
  1871.             else
  1872.             {
  1873.               return(TRUE);
  1874.             }
  1875.           }
  1876.         }
  1877.  
  1878.         error--;
  1879.       }
  1880.       else
  1881.       {
  1882.         if (error > -EPSILON)
  1883.         {
  1884.           z += signz;
  1885.  
  1886.           error--;
  1887.         }
  1888.       }
  1889.  
  1890.       x += signx;
  1891.  
  1892.       error += delta;
  1893.  
  1894.       neary = fary;
  1895.  
  1896.       fary += deltay;
  1897.     }
  1898.     while ((neary >= ymin) && (neary <= ymax) && (x >= xmin) && (x <= xmax) && (z >= zmin) && (z <= zmax));
  1899.   }
  1900.   else
  1901.   {
  1902.     /*
  1903.      * Z-axis is driving axis.
  1904.      */
  1905.  
  1906.     delta = fabs(dx / dz);
  1907.  
  1908.     x = (int)px;
  1909.     z = (int)pz;
  1910.  
  1911.     x0 = px - floor(px);
  1912.     z0 = pz - floor(pz);
  1913.  
  1914.     signx = sign(dx);
  1915.     signz = sign(dz);
  1916.  
  1917.     /* Get initial error. */
  1918.  
  1919.     if (dz >= 0.0)
  1920.     {
  1921.       if (dx >= 0.0)
  1922.       {
  1923.         error = x0 + delta * (1.0 - z0) - 1.0;
  1924.       }
  1925.       else
  1926.       {
  1927.         error = -(x0 - delta * (1.0 - z0));
  1928.       }
  1929.     }
  1930.     else
  1931.     {
  1932.       if (dx >= 0.0)
  1933.       {
  1934.         error = x0 + delta * z0 - 1.0;
  1935.       }
  1936.       else
  1937.       {
  1938.         error = -(x0 - delta * z0);
  1939.       }
  1940.     }
  1941.  
  1942.     /* Get y differential. */
  1943.  
  1944.     deltay = dy / fabs(dz);
  1945.  
  1946.     if (dz >= 0.0)
  1947.     {
  1948.       neary = Start[Y] - z0 * deltay;
  1949.  
  1950.       fary = neary + deltay;
  1951.     }
  1952.     else
  1953.     {
  1954.       neary = Start[Y] - (1.0 - z0) * deltay;
  1955.  
  1956.       fary = neary + deltay;
  1957.     }
  1958.  
  1959.     /* Step through the cells. */
  1960.  
  1961.     do
  1962.     {
  1963.       if (neary < fary)
  1964.       {
  1965.         y1 = neary;
  1966.         y2 = fary;
  1967.       }
  1968.       else
  1969.       {
  1970.         y1 = fary;
  1971.         y2 = neary;
  1972.       }
  1973.  
  1974.       if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1975.       {
  1976.         if (HField->Type & IS_CHILD_OBJECT)
  1977.         {
  1978.           found = TRUE;
  1979.         }
  1980.         else
  1981.         {
  1982.           return(TRUE);
  1983.         }
  1984.       }
  1985.  
  1986.       if (error > EPSILON)
  1987.       {
  1988.         x += signx;
  1989.  
  1990.         if ((x < xmin) || (x > xmax))
  1991.         {
  1992.           break;
  1993.         }
  1994.         else
  1995.         {
  1996.           if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1997.           {
  1998.             if (HField->Type & IS_CHILD_OBJECT)
  1999.             {
  2000.               found = TRUE;
  2001.             }
  2002.             else
  2003.             {
  2004.               return(TRUE);
  2005.             }
  2006.           }
  2007.         }
  2008.  
  2009.         error--;
  2010.       }
  2011.       else
  2012.       {
  2013.         if (error > -EPSILON)
  2014.         {
  2015.           x += signx;
  2016.  
  2017.           error--;
  2018.         }
  2019.       }
  2020.  
  2021.       z += signz;
  2022.  
  2023.       error += delta;
  2024.  
  2025.       neary = fary;
  2026.  
  2027.       fary += deltay;
  2028.     }
  2029.     while ((neary >= ymin-EPSILON) && (neary <= ymax+EPSILON) &&
  2030.            (x >= xmin) && (x <= xmax) &&
  2031.            (z >= zmin) && (z <= zmax));
  2032.   }
  2033.  
  2034.   return(found);
  2035. }
  2036.  
  2037.  
  2038.  
  2039. /*****************************************************************************
  2040. *
  2041. * FUNCTION
  2042. *
  2043. *   block_traversal
  2044. *
  2045. * INPUT
  2046. *
  2047. *   Ray    - Current ray
  2048. *   HField - Height field
  2049. *   Start  - Start point for the walk
  2050. *
  2051. * OUTPUT
  2052. *
  2053. * RETURNS
  2054. *
  2055. *   int - TRUE if intersection was found
  2056. *   
  2057. * AUTHOR
  2058. *
  2059. *   Dieter Bayer
  2060. *
  2061. * DESCRIPTION
  2062. *
  2063. *   Traverse the blocks of the height field.
  2064. *
  2065. * CHANGES
  2066. *
  2067. *   Feb 1995 : Creation.
  2068. *
  2069. *   Aug 1996 : Fixed bug as reported by Dean M. Phillips:
  2070. *              "I found a bug in the height field code which resulted
  2071. *              in "Illegal grid value in dda_traversal()." messages
  2072. *              along with dark vertical lines in the height fields.
  2073. *              This turned out to be caused by overlooking the units in
  2074. *              which some boundary tests in two different places were
  2075. *              made. It was easy to fix.
  2076. *
  2077. ******************************************************************************/
  2078.  
  2079. static int block_traversal(Ray, HField, Start)
  2080. RAY *Ray;
  2081. HFIELD *HField;
  2082. VECTOR Start;
  2083. {
  2084.   int xmax, zmax;
  2085.   int x, z, nx, nz, signx, signz;
  2086.   int found = FALSE;
  2087.   int dx_zero, dz_zero;
  2088.   DBL px, pz, dx, dy, dz;
  2089.   DBL maxdv;
  2090.   DBL ymin, ymax, y1, y2;
  2091.   DBL neary, fary;
  2092.   DBL k1, k2, dist;
  2093.   VECTOR nearP, farP;
  2094.   HFIELD_BLOCK *Block;
  2095.  
  2096.   px = Start[X];
  2097.   pz = Start[Z];
  2098.  
  2099.   dx = Ray->Direction[X];
  2100.   dy = Ray->Direction[Y];
  2101.   dz = Ray->Direction[Z];
  2102.  
  2103.   maxdv = (dx > dz) ? dx : dz;
  2104.  
  2105.   /* First test for 'perpendicular' rays. */
  2106.  
  2107.   if ((fabs(dx) < EPSILON) && (fabs(dz) < EPSILON))
  2108.   {
  2109.     x = (int)px;
  2110.     z = (int)pz;
  2111.  
  2112.     neary = Start[Y];
  2113.  
  2114.     if (dy >= 0.0)
  2115.     {
  2116.       fary = 65536.0;
  2117.     }
  2118.     else
  2119.     {
  2120.       fary = 0.0;
  2121.     }
  2122.  
  2123.     return(intersect_pixel(x, z, Ray, HField, min(neary, fary), max(neary, fary)));
  2124.   }
  2125.  
  2126.   /* If we don't have blocks we just step through the grid. */
  2127.  
  2128.   if ((HField->Data->block_max_x <= 1) && (HField->Data->block_max_z <= 1))
  2129.   {
  2130.     return(dda_traversal(Ray, HField, Start, &HField->Data->Block[0][0]));
  2131.   }
  2132.  
  2133.   /* Get dimensions of grid. */
  2134.  
  2135.   xmax = HField->Data->block_max_x;
  2136.   zmax = HField->Data->block_max_z;
  2137.  
  2138.   ymin = (DBL)HField->Data->min_y - EPSILON;
  2139.   ymax = (DBL)HField->Data->max_y + EPSILON;
  2140.  
  2141.   dx_zero = (fabs(dx) < EPSILON);
  2142.   dz_zero = (fabs(dz) < EPSILON);
  2143.  
  2144.   signx = sign(dx);
  2145.   signz = sign(dz);
  2146.  
  2147.   /* Walk on the block grid. */
  2148.  
  2149.   px /= HField->Data->block_width_x;
  2150.   pz /= HField->Data->block_width_z;
  2151.  
  2152.   x = (int)px;
  2153.   z = (int)pz;
  2154.  
  2155.   Assign_Vector(nearP, Start);
  2156.  
  2157.   neary = Start[Y];
  2158.  
  2159.   /*
  2160.    * Here comes the block walk algorithm.
  2161.    */
  2162.  
  2163.   do
  2164.   {
  2165. #ifdef HFIELD_EXTRA_STATS
  2166.     Increase_Counter(stats[Ray_HField_Block_Tests]);
  2167. #endif
  2168.  
  2169.     /* Get current block. */
  2170.  
  2171.     Block = &HField->Data->Block[z][x];
  2172.  
  2173.     /* Intersect ray with bounding planes. */
  2174.  
  2175.     if (dx_zero)
  2176.     {
  2177.       k1 = BOUND_HUGE;
  2178.     }
  2179.     else
  2180.     {
  2181.       if (signx >= 0)
  2182.       {
  2183.         k1 = ((DBL)Block->xmax + 1.0 - Ray->Initial[X]) / dx;
  2184.       }
  2185.       else
  2186.       {
  2187.         k1 = ((DBL)Block->xmin - Ray->Initial[X]) / dx;
  2188.       }
  2189.     }
  2190.  
  2191.     if (dz_zero)
  2192.     {
  2193.       k2 = BOUND_HUGE;
  2194.     }
  2195.     else
  2196.     {
  2197.       if (signz >= 0)
  2198.       {
  2199.         k2 = ((DBL)Block->zmax + 1.0 - Ray->Initial[Z]) / dz;
  2200.       }
  2201.       else
  2202.       {
  2203.         k2 = ((DBL)Block->zmin - Ray->Initial[Z]) / dz;
  2204.       }
  2205.     }
  2206.  
  2207.     /* Figure out the indices of the next block. */
  2208.  
  2209.     if ((k1 < k2 - EPSILON / maxdv) && (k1 > 0.0))
  2210.     {
  2211.       /* Step along the x-axis. */
  2212.  
  2213.       dist = k1;
  2214.  
  2215.       nx = x + signx;
  2216.       nz = z;
  2217.     }
  2218.     else
  2219.     {
  2220.       if ((k1 < k2 + EPSILON / maxdv) && (k1 > 0.0))
  2221.       {
  2222.         /* Step along both axis (very rare case). */
  2223.  
  2224.         dist = k1;
  2225.  
  2226.         nx = x + signx;
  2227.         nz = z + signz;
  2228.       }
  2229.       else
  2230.       {
  2231.         /* Step along the z-axis. */
  2232.  
  2233.         dist = k2;
  2234.  
  2235.         nx = x;
  2236.         nz = z + signz;
  2237.       }
  2238.     }
  2239.  
  2240.     /* Get point where ray leaves current block. */
  2241.  
  2242.     VEvaluateRay(farP, Ray->Initial, dist, Ray->Direction);
  2243.  
  2244.     fary = farP[Y];
  2245.  
  2246.     if (neary < fary)
  2247.     {
  2248.       y1 = neary;
  2249.       y2 = fary;
  2250.     }
  2251.     else
  2252.     {
  2253.       y1 = fary;
  2254.       y2 = neary;
  2255.     }
  2256.  
  2257.     /* Can we hit current block at all? */
  2258.  
  2259.     if ((y1 <= (DBL)Block->ymax + EPSILON) && (y2 >= (DBL)Block->ymin - EPSILON))
  2260.     {
  2261.       /* Test current block. */
  2262.  
  2263. #ifdef HFIELD_EXTRA_STATS
  2264.       Increase_Counter(stats[Ray_HField_Block_Tests_Succeeded]);
  2265. #endif
  2266.  
  2267.       if (dda_traversal(Ray, HField, nearP, &HField->Data->Block[z][x]))
  2268.       {
  2269.         if (HField->Type & IS_CHILD_OBJECT)
  2270.         {
  2271.           found = TRUE;
  2272.         }
  2273.         else
  2274.         {
  2275.           return(TRUE);
  2276.         }
  2277.       }
  2278.     }
  2279.  
  2280.     /* Step to next block. */
  2281.  
  2282.     x = nx;
  2283.     z = nz;
  2284.  
  2285.     Assign_Vector(nearP, farP);
  2286.  
  2287.     neary = fary;
  2288.   }
  2289.   while ((x >= 0) && (x < xmax) && (z >= 0) && (z < zmax) && (neary >= ymin) && (neary <= ymax));
  2290.  
  2291.   return(found);
  2292. }
  2293.  
  2294.